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ABSTRACT 


An  algorithm  for  the  identification  of  the  input  distribution 
matrix  of  a  linear3  stationary  system  operating  in  a  stochastic 
environment  is  derived.   The  identification  is  accomplished  by  defining 
a  set  of  autocorrelation  functions  for  a  noise  element  composed  of  a 
linear  combination  of  the  input  distribution  matrix  elements  and  the 
random  excitations  of  the  system.   Another  possible  identification 
method  employing  a  Kalman  filter  is  discussed  and  the  problems  associ- 
ated with  its  derivation  are  presented.   Results  of  computer  simulations 
for  both  methods  are  included. 
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I.   INTRODUCTION 

In  recent  years,  an  emphasis  on  plant  identification  for  a  system 
operating  in  a  stochastic  environment  has  developed  in  the  field  of 
automatic  control.   It  is  necessary  to  know  the  dynamics  of  the  system 
before  it  may  be  effectively  controlled.   In  order  to  know  the  system 
dynamics,  a  suitable  mathematical  model  of  the  system  must  be  found. 
For  a  system  subjected  to  a  random  forcing  function,  finding  a  suitable 
mathematical  model  involves  the  identification  of  both  the  state 
transition  matrix  and  the  input  distribution  matrix  for  the  random 
forcing  function.   Several  methods  have  been  proposed  for  the  identi- 
fication of  the  state  transition  matrix,  but  no  one  has  proposed  a 
method  for  identifying  the  input  distribution  matrix  for  systems'  with 
a  single  random  input.   The  purpose  of  this  investigation  is  to  propose 
such  a  method.   The  problem  may  be  defined  as  follows: 

Given : 

1.  A  system  whose  behavior  may  be  defined  by  a  set  of  linear, 
constant-coefficient,  differential  equations. 

2.  A  statistical  description  of  the  excitation  of  the  system. 

3.  A  sequence  of  observations  on  the  state  vector. 

4.  An  estimate  of  the  state  transition  matrix  of  the  system. 
Problem: 

Determine  the  input  distribution  matrix  of  the  system. 
Although  this  investigation  is  concerned  only  with  systems  with 
random  forcing  functions,  the  methods  described  may  be  extended  to 
systems  with  both  random  and  control  inputs.   The  method  of  section 
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Ill  may  be  extended  by  changing  the  form  of  the  estimated  state  transi- 
tion matrix  as  described  by  Lee  [5];  and  that  of  section  IV,  by  including 
the  effect  of  the  control  input  in  the  filter  equations  as  described  by 
Kalman  [3].   The  problem  is  restricted  to  random  forcing  functions  for 
simplicity. 
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II.   BASIC  MATHEMATICAL  MODEL 

The    investigation   reported   here    starts   with    the   assumption   that 
there    exists    a    set    of    n    linear,    constant-coefficient,    differential 
equations    of   the    form 

x(t)  -  Ax(t)  +  Bw(t)  (2.1) 

which  describe  the  dynamic  behavior  of  a  system  subjected  to  a  random 
forcing  function,  w(t).   The  solution  of  these  equations  is 

x(t)  =  $(t-tQ)  x(0)  +  ]    $(t-T)  Bw(T)dT  (2.2) 

where 


Because  the  continuous  system  dynamics  will  be  observed  at  discrete 
times,  it  is  necessary  to  describe  the  system  by  discrete  or  difference 
equations.   Introducing  a  sampling  device  of  period  T  and  a  zero- 
order  hold  on  the  forcing  function  allows  the  system  equations  to  be 
written 

2£k+i  =  $(T)2£k  +  r('r)^k>  <2'3> 

where   Xk    ^s    tne   system   output   vector   at   a    time   kT,    $    is    the    nxn  discrete 

state  transition  matrix,  T(T)  is  the  lxn  input  distribution  matrix  and 

w      is    the    sampled   and    zero-order    held    input   vector   at    time   kT.      The 

components  of  the  input  vector  w,  are  assumed  to  form  a  gauss ian  white 

sequence  of  zero  mean  and  known  variance  during  the  period  of  system 

operation  under  investigation. 

Observations  made  on  the  system  output  are  assumed  to  be  scalar 

linear  combinations  of  the  states  of  the  system  which  are  contaminated 

by  additive  gauss ian  white  noise  of  zero  mean  and  known  variance.   At 
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the   k        sampling    instant    the    observation   z      is   given  as 

K. 

zk  =  H^k  +  vk 

where  H  is  the  lxn  known  measurement  matrix  and  v,  is  the  scalar 

k 


(2.4) 


additive  measurement  noise.   The  system  may  then  be  represented  by 
the  block  diagram  of  Figure  1. 

In  the  general  problem  of  identification  of  linear  dynamic  systems, 
it  is  necessary  to  estimate  the  elements  of  the  state  transition  matrix 
before  proceeding  to  identify  the  input  distribution  matrix.   For  the 
purposes  of  this  investigation,  this  estimation  is  assumed  done  by 
using  the  method  described  by  Lee  [5],   The  estimation  of  the  state 
transition  matrix  by  this  method  makes  it  necessary  to  change  the 
system  model  to  the  so-called  canonic  form.   Assuming  that  the  system 
is  observable  and  n-ident if iable  ,  the  difference  equation  for  the 
system  model  then  becomes 


^+1  -  **  ^  +  r*  ^ 


\  =  H*  ^k  +  \ 


(2.5) 


where 


-k+1 


C   = 


=   C 


\+l 


H 

H$  . 


H$ 


n-1 


$*  =  C$C 


-1 


0  ! 
?       i 

.     I 

•     I 

Cp1    CD2     .  .  . 

cp 
n 

(2.6) 


14 


r*  =  cr  =  b  = 


H*  =  [  1  0  0  •  •  •  0  ] 

The  matrix  is  now  the  estimated  state  transition  matrix  obtained  by 
using  the  method  of  Lee.   The  canonic  system  model  can  then  be 
represented  by  the  block  diagram  of  Figure  2„ 
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FIGURE  1.   BLOCK  DIAGRAM  OF  THE  SYSTEM 


w. 


r* 


FIGURE  2.   BLOCK  DIAGRAM  OF  THE 

SYSTEM  IN  CANONIC  FORM 
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III.   DERIVATION  OF  THE  IDENTIFICATION  ALGORITHM 
Consider  a  system  with  canonic  state  equations  of  the  form 

ik+l  =  **  ik  +  r*  wk 

(3.1) 

Zk+1  =  H*  *k+l  +  \ 

where  w,  is  assumed  to  be  a  scalar  random  forcing  function.   The 
— k 

z-trans formed  state  equation  becomes 

z£(z)    -   z£(0)   =  $*  ^(z)  +  T*  w(z).  (3.2) 

Assuming    that    the    initial   conditions  _y_(0)    are   zero   or    that    the 
identification   process    does    not    start    until    the    initial   conditions 
no    longer   affect    the    system  output,    the    system  may   be    represented   by 
the    signal-flow  graph   of   Figure    3.      From   this    signal-flow  graph,    the 
transfer    function  Z{z)/w(z)   may   be   written  as 

b  b       .    ^  Cpxb  cp  <P      -r\ 

„ ,   v  n  n-1  v  z  y  n-2   v         z  2  y 

Z(z)   _  z  z z z  -3   3> 

W(z)  cp         cp      ,  cd0  cp, 

1        _J1  n- 1        1 1  t  2  _1 

z  2  n-1  n 

z  z  z 

Since  the  elements  of  the  $*  matrix  are  assumed  to  be  known,  the 

problem  can  now  be  stated  as  finding  the  location  of  the  zeros  of 

the  transfer  function  and  thus  the  values  of  b, ,  b„ ,  ...,  b  . 

12  n 

From   the    transfer    function,    the   difference    equation   relating    the 
measured    output    z      to    the    input   w     may   be   written  as 

K.  K. 
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z.     =  v  +  cp  z,       +  cp0z         ,!+-..+  Cp  z,  -b.w.       -(b0-cp  b^w,     _ 
k  k      M   k-n     ^2   k-n+1  ^n  k      1   k-1        2   ^n    1      k-2 


(b3-cpnb2-cp^1b1)wk_3 


(3.4) 


-...-(b      ,-cpb      _-...-co_b    )w 

n-1      n   n-2  3    1     k-n+1 


(b    -cp  b      ,-.  .  .-cp  b    )w 
n^nn-1  2    1     k-n 


Shifting  indices  by  +1  and  using  vector  notation,  Equation  (3.4) 

T 
becomes  zfc+1  =J     eg  +  n^d  +  vR  (3.5) 

—  k 


where 


1- 


'k-n+1 
:k-n+2 


L>    J 


C£  = 


CPr 


^ 


d   = 


wk 

V 

•  1 

V 

-n+1 

— 

_ 

'-I  x 


\ 


-cp 

"CP. 


1   X 


-cp 
n- 1     ^n 


-op. 


-<P- 


\ 


1\ 


»  \ 


-cp. 


(3.6) 


r-       -* 

bl 

= 

b2 

\ 

* 

1 
__ 

b 
n 

^_   _ 
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Because  the  elements  of  both  the  vectors  n  and  d  are  unknown,  a  single 

~~k  —  ° 

noise    element    is    defined    as 

Ok   =  nkTd  +  vk  (3.7) 


where   Q,  is    a    scalar   and 

Q.  =  d,    w.    +  d„   w,     .+....+  d      w,        .  ..   +  v, 
k  Ik  2      k-1  n     k-n+1  k 

Q     ,  =  d,    w  +  d0   w     +    . . .   +  d     w,  0  +  v,  ,  , 

Tc+1  1      k+1  2      k  n     k-n+2  k+1 


(3.8) 


<W   =  dl  \+n  +  d2   Viri-1  +    •  •'   +  dn  Wk+1  +  \+n 

From  Equation    (3.8)    it    is    apparent    that    the    noise   element  Q     is 

correlated   with    the   elements    Q     .,    ^Vio  >    •••»    £\i       ■,  ;    however,    0 

and   Q         are    independent  with   respect    to   each   other.      This    property 

of   the    noise   elements    is    the   basis    for   the    identification   algorithm, 

As    shown   in  Appendix  A,    a    set    of  autocorrelation   functions    for    the 

noise   elements   may   be   defined    such    that 

2,  2    ,,   2         .   2  2S  2 


R 


(0)  =E   [^  ]  =aw     (dL    +d2    +  ...  +dn  )  +  a^ 


(3.9) 


2   n"- 


-   1,    2,    . . . ,    n-1 


R(j)  =  E   [cy^]  =  aw     £  d.  dk+.  j  -  1 

R(j)    =0  j   >  n 

where  a   is  the  standard  deviation  of  the  gaussian  white  sequence  of 
w 

noise  inputs   w,  and  a   is  the  standard  deviation  of  the  measurement 
r      k      v 


noise . 


By   rearranging   Equation    (3.5)    such    that 

Z  k 


T 

\   =  \+l   ~%.    eg    ,  (3.10) 
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the    noise    element    is    now  a    linear   function   of   the    system   outputs 

z,  ,  ,    z,        ,~,    .  .  .  ,    z,  ,  ,    and    is    therefore   known.      Hence,    the   auto- 

k-n+1        k-n+2  k+1 

correlation   functions    of  equation    (3.9)   may   be    found   by    taking    the 

2 
statistics    E    [fl     ],    E    [Qfl      .].      Equation    (3.9)    now  represents   a    set 

of   n   simultaneous    nonlinear   algebraic   equations   which  may   be    solved 

for    the    values    of  d, ,    d^ ,     ....    d    .      The   values    of  b,,    b„ ,    ....    b     may 

12  n  12  n 

then  be    found    from  equation    (3.6). 

The   general   procedure    for   the    identification   of   the    input    distri- 
bution matrix   in   a    linear   system  with   a    random   forcing    function  may  be 
outlined   as    follows: 

1.  Measure    a   scalar    linear   combination   of   the    states    of    the    system 
plus    additive    noise. 

2.  Identify   the   canonic    form  of    the    state    transition  matrix    of 
the    system. 

3.  Calculate    the   values    of   the    noise    element    fi     at   each   sampling 

K. 

time . 

4.  Calculate  the  autocorrelation  functions  R(j)  for  j  =0,  1,  ..., 
n-1. 

5.  Solve  the  simultaneous  equations  (3.9)  for  the  elements  of  d. 

6.  Solve  equation  (3.6)  for  the  elements  of  b_. 

Although  the  above  algorithm  does  provide  a  means  for  identifying 
the  input  distribution  matrix,  two  problems  arise  in  the  identification 
process.   The  first  of  these  involves  deciding  when  to  start  the 

identification  process.   Because  the  autocorrelation  functions  are 

2 
computed  by  taking  the  average  of  the  values  of  Q  ,  0,0,,-,,  ..., 

QQ,    ..  at  a  number  of  sampling  instants,  this  process  must  not  begin 

until  enough  time  has  elapsed  so  that  the  sample  mean  and  variance  of 
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2  2 

the  noise  input  and  the  measurement  noise  approach  0,  a   and  0,  a 

w  v 

respectively.   If  the  identification  process  is  begun  too  soon,  the 
computed  values  of  _b  will  be  in  error  from  the  actual  values. 

The  other  problem  in  the  identification  process  can  be  best  shown 
by  an  example.   Assuming  n  =  2,  the  autocorrelation  functions  then 

be  c  ome 

2    2     2      2 
R(O)  -a    (d/  +  d0  )  +  o 

w     L      I  v 


(3.11) 


K(D  -  aw  dl  d2  . 

Solution   of   Equation    (3.11)    gives 


dl=±L 


1 


2a 


2  - 


2,2 


2  |_R(0)    -  av     +(JR(0)    -  cv'r   -  4[R(1)£] 


l/2n 


1/2 


w 


(3.12) 


±B 


2       —  L        2 
L2a 
w 


R(0)   -  o^2  ±  Q[R(o)   -  av2]2 


-   4[R(l)2f) 


1/2-,    -n  1/2 


and  from  Equation  (3.6) 


bl^dl 


b2   "  d2  +  cP2bl 


(3.13) 


It    is    apparent    from  Equation    (3.12)    that    the    identification   of    the 
elements    b, ,    b„    involves    a   decision  as    to  which    combination   of   signs 
to   use    as    the    correct    one.      For   higher    order   systems    this   decision 
becomes    more    complex   since    the    number    of    possible    sign   combinations 
increases . 
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IV.   KALMAN  FILTER  APPROACH 

An  alternative  to  the  identification  algorithm  of  the  preceding 
section  may  be  derived  from  the  equations  of  the  Kalman  filter 
described  below.   Consider  again  a  system  with  the  canonic  state 
equations  of  the  form 

*k+l  =  **\   +  r*  ^k 

(4.1) 

zk+l   =  H*  *k+l  +  Vl 
where   w,     is    a    scalar   random   forcing    function.      The    canonic    state 
vector   v,     may   be    estimated   by   the    technique    described   by  Kalman    [3] 

K. 

where,  given  k  previous  measurements,  the  estimate  jv /,  at  the  k 
sampling  instant  is  given  by 

4/k =  **  i-i/k-i +  Gk  [zk  - H*  **  4-i/k-i1         <4-2> 

where    G,     is    the    filter  weighting    or   gain  applied   at    the   k        iteration 
and   H*   is    the    lxn  known  measurement   matrix.      The   gain    is    calculated 
from  a    set    of    recursive    equations    of   the    form 

Gk   "    Pk/k-l   H+T    IH*   Vk-1   H*T  +  Rl"1 

Pk/k   "    (1   '   V*    Pk/k-l  (4.3) 

Pk+l/k   "   **  Pk/k   **T  +  Q 
where    P    .      and    P         ,      are    the    estimation  and    prediction  covariances, 
respectively,    and   R  and  Q   are   given  by 


A  2 

R   =  E    [vRZ] 

A  2  T 

q  =  r*  e  [w.    ]  r*1 


(4.4) 
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Since  all  elements  necessary  to  predict  the  canonic  state  vector 
v,  are  known  except  the  input  distribution  matrix  T*,  the  system  model 
based  on  this  estimation  scheme  will  produce  suitable  estimates  for  y_ 

K. 

only  if  the  correct  T*   is  used  in  the  filter  equations.   In  order  to 
determine  if  the  correct  F*  is  being  used,  it  is  necessary  to  define  a 
performance  index  to  compare  the  performance  of  the  model  with  that  of 
the  system.   Because  the  actual  system  output  is  available  only  in  the 
form  of  the  measurement  z.  ,  the  performance  index  must  be  defined  in 
terms  of  z   and  the  predicted  measurement  z  ,      ,  as  given  by 

Vk-irBM*ti/k-i'  (4-5) 

which    is    obtained    from   the    system  model. 

Before   defining    the    performance    index   to   be    used,    it    is    necessary 
to   define    some    other    terms.      The    filter    is    said   to   have    reached    the 
steady-state    condition  when   the    prediction  covariance   matrix  does    not 
change    from   one    sampling    instant    to    the    next.      Under    the    steady-state 
condition   the   actual    prediction  measurement  variance    is   given   by 

,    K+M 
EZ(M)    --     S      (z.    -    z\  )  (4.6) 

i  =  K 

where  K  is  the  value  of  k  at  the  sampling  instant  when  the  filter 

reaches  the  steady-state  condition.   The  integer  M  is  chosen  such 

that  the  relation 

|EZ(M+1)  -  EZ(M)|  <  0.05|EZ(M)|  (4.7) 

is  valid.   The  predicted  measurement  variance  is  defined  as 

PZ  =  EK\  -  Vk-i)2'  <4-8> 


24 


which,  as  shown  in  Appendix  B,  may  be  written  as 


PZ  =  H*  P H*  +  R 
K 


(4.9) 


where    P      is    the    prediction   covariance    under    the    steady-state    conditions 
K 

of   the    filter. 

From   the   above   definitions    the    performance    index  J   is    defined   as 
the   difference    between   the   actual   and    the    predicted   measurement 
variances    and   may   be   written  as 

J   =  EZ(M)    -    PZ 


J   = 


-,    K+M 
-       L        (Z.     -    Z.  ,  ) 

LM    .   „        1  l/i-l 

l-K 


H*   P^   H*     +   R 
K 


(4.10) 


Since  it  is  desired  that  the  output  of  the  model  z  .    and  the  output 

of  the  system  z   be  equal  for  perfect  operation  of  the  model,  J  is  to 

be  minimized.   Because  both  EZ(M)  and  PZ  are  functions  of  the  input 

distribution  matrix  T*,  the  minimum  for  J  will  result  when  the  correct 

r*  is  used  in  the  model.   The  problem  may  now  be  stated  for  this 

approach  as  identifying  F*  so  as  to  minimize  the  performance  index  J. 

Because  it  would  waste  time  and  effort,  a  tria 1-and-error  search 

for  the  correct  V*   is  undesirable.   Therefore,  it  is  necessary  to 

derive  some  method  of  changing  the  elements  of  T*  used  in  the  model 

after  each  calculation  of  J,  based  on  the  value  of  J.   The  method 

investigated  is  based  on  a  gradient  search  where  the  elements  of  T* 

are  changed  by  some  amount  according  to  the  values  of  a  sensitivity 

function  and  J.   The  sensitivity  function  must  reflect  the  change  in 

J  caused  by  changing  T*.   The  method  of  changing  T*  may  be  written  as 

r*   =  T*   +  f  (J,  r*)  J  (4.11) 

new    old   — 
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where  the  subscripts  "old"  and  "new"  represent  the  value  of  T*  just 

used  in  the  model  and  the  value  to  be  tried  next,  respectively.   The 

sensitivity  function  is  denoted  by  _f  (J,  T*)  where  f  is  a  vector 

function. 

The  problem  is  now  to  find  a  suitable  sensitivity  function  to  use 

Two  functions  and  their  failures  during  simulation  are  presented  in 

the  next  section.   Because  the  performance  index  is  a  function  of  the 

measurements  z  „  z„ „ n  ,  ....  z„lw  (which  are  functions  of  all  past 
K   Kfl        K+M  r 

measurements)  and  the  predicted  measurements  during  the  interval  K, 
K+M,  (which  are  also  functions  of  all  past  measurements),  the  problem 
of  finding  a  suitable  sensitivity  function  is  a  complicated  one.   At 
present  no  suitable  function  has  been  found. 
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V.   COMPUTER  SIMULATIONS 


In  order  to  determine  the  accuracy  of  the  proposed  identification 
algorithm,  a  second-order  system  was  simulated  on  the  IBM  360/67  com- 
puter.  For  the  system  chosen  the  canonic  state  transition  matrix,  the 
input  distribution  matrix,  and  the  measurement  matrix  were 

0      1 
-1/16   1/2 

1 
1/2 


i*  = 


r*  = 


(5.1) 


H* 


■[»"]. 


In   the    simulation,    the   actual   and    the    canonic-state    transition  matrices 
were   assumed   equal    so    that    no   error    in   the    identification   process   would 
result    from  a   difference    in   the    two^and   also    for    simplicity.      The    random 
forcing    function  and    the   measurement    noise   were    obtained    from  a    random 
number  generator  which  gives    a   gaussian  white    sequence    of  any   desired 
mean  and    standard   deviation.      The   means    and    standard   deviations    for   the 


two   sequences   were 


E    [wk]    =0 
E    [wk2 ]    =  2 
E    [vk]    =   0 
E    [vk2]    .    1 


The    system  dynamics   were    simulated   and    the   measurements    z,     were    used 

k 

to   obtain   the   autocorrelation   functions    R(0)    and   R(l).      These   values 


(5.2) 
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were  'Chen   used    to   compute    the   values    of   b1    and   b_    at   each   sampling 
instant    from    the    equations 

1/2-,  -il/2 


^k)    =[ 


0.125 


R(0)    -    1  +    [(R(0)    -    l)2    -   4R(1)2] 


(5.3) 


b2(k)    = 


0.125 


R(0)    -    1    -    [(R(0)    -    l)2    -   4R(1)2] 


1/2-,  I 


1/2 


The    system  and   the    identification   process   were   simulated    for    1000 
sampling    intervals    of    period    one    second.      An   ensemble    of   200   simu- 
lations  was    run  and   an   ensemble   average    for   R(0) ,    R(l)  ,    b..(k),    and 
b~(k)   was    taken.      The    results    at   various   sampling    intervals    are    shown 
in  Table    I.      The   ensemble    averages    of    the   differences    between    the 
actual   values    of  F*  and    the    computed   values   are   also   shown. 

As    can  be    seen   from  Table    I,    there    is    a   bias    in   the   autocorrelation 
function  R(l)   which    produces    a   bias    in   the   estimates    for   b      and   b„ . 
No   explanation    for    the    bias    can   be   given.      Simulation    of  another, 
similar   second-order   system  also   produced   a   bias,    equal    to    -1.0. 
For    the   algorithm   to    produce    the   desired   results,    it    is    necessary   to 
eliminate    the    bias. 

A   second-order   system  was    also   simulated    for    the   Ka lman-f ilter 
approach    to    the    identification   problem.      The   values    used    for   the 
simulation  were 

0  1 

1/8      1/4| 
1 


>*  = 


r*  = 

1/2 

—     

H*  =     fl      0 


(5.4) 
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The    random   forcing    function  and    the   measurement    noise  were   generated   as 
described   above   with 


E    [wfc]    =   0 

E    [wk2]    =   3 

E    [vk]    =   0 

E    [v,  2]    =  R  =  4     . 
k 


(5.5) 


The  performance  index  J  was  calculated  for  various  values  of  T*   with 
the  results  given  in  Table  II.   It  can  be  seen  that  J  is  a  minimum 
when  the  correct  T*   is  used. 

Two  sensitivity  functions  were  also  calculated  during  the  simulation 
The  first  as  derived  in  Appendix  C  is  given  by 


k  v-  r*>  ■  ft 


The  second,  also  derived  in  Appendix  C,  is  given  by 


(5.6) 


i2  (J,  r*>  =  h*  ^ 


dG   r 

K 


a      21    ^PZ 

zv  (zi  "  zi/i-i}  J  "  W* 


(5.7) 


As  shown  in  Table  II,  neither  of  these  sensitivity  functions  is 
suitable  since  neither  produces  a  consistent  pattern  by  which  to 
change  r*.   The  failure  of  these  functions  appears  to  lie  in  the  fact 
that  neither  measures  the  sensitivity  of  the  measurements  to  the  values 

of  r*. 
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bL(k)  1-b^k)  b2(k)  l/2-b2(k)        E[^2]k   E[nkQk+1]k 


2  0.0  0.0  0.0  0.0  5.11  0.0 

10  0.84  0.16  0.65  -0.15  5.25  -0.31 

100  1.003  -0.003  0.645  -0.145  5.17  -0.49 

200  1.009  -0.009  0.635  0.135  5.18  -0.50 

300  1.010  -0.01  0.633  -0.133  5.17  -0.50 

400  1.012  -0.013  0.633  -0.133  5.18  -0.50 

500  1.016  -0.016  0.631  -0.131  5.20  -0.50 

600  1.017  -0.017  0.634  -0.134  5.21  -0.51 

700  1.016  -0.016  0.633  -0.133  5.20  -0.51 

800  1.018  -0.018  0.634  -0.135  5.22  -0.51 

900  1.019  -0.019  0.636  -0.136  5.23  -0.51 

999  1.019  -0.019  0.635  -0.135  5.23  -0.51 


2  2  2  2  2 

R(O)    =  a        Cd-,      +  d„    )   +  a        =   5.0    (should  match   E[Q     ]) 
w      v    1  2  v  k 

R(l)    =  aw     d   d2    =  0.0    (should   match   E^Q^^]) 


TABLE   I.      Results    of   Simulation 
of   Identification  Algorithm 
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hv 

b2 

EZ(M) 

PZ 

J 

iL  (J  = 

r*) 

i2  (J  = 

,  r*) 

1/2, 

1/4 

14.7 

6.7 

8.7 

-8.4, 

-2.3 

-2.4, 

-0.6 

1, 

1/4 

14.6 

13.4 

1.2 

-17.7, 

-0.3 

12.1, 

0.2 

3/2, 

1/4 

15.0 

24.8 

-9.8 

-26.6, 

0.4 

24.1, 

-0.3 

1/2, 

1/2 

14.1 

7.4 

6.7 

-6.3, 

-4.2 

0.1, 

0.04 

*1, 

1/2 

13.7 

13.3 

0.4 

-18.0, 

1.1 

12.5, 

-0.8 

3/2, 

1/2 

14.3 

24.5 

-11.2 

-27.3, 

2.1 

24.7, 

-1.9 

1/2, 

3/4 

15.7 

8.5 

7.2 

-0.3, 

-7.6 

0.1, 

1.3 

1, 

3/4 

19.9 

11.0 

8.9 

-60.5, 

60.6 

23.0,- 

■23.0 

3/2, 

3/4 

15.0 

23.0 

-8.0 

-31.6, 

8.9 

28.1, 

-7.9 

*Correct  T* 


TABLE  II.   Results  of  Simulation  of 
Kalman  Filter  Approach 
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VI.   CONCLUSIONS 

The  identification  of  the  input  distribution  matrix  of  a  linear 
system  operating  in  a  stochastic  environment  is  a  complicated  process, 
but  it  can  be  done  to  some  degree  by  the  proposed  algorithm.   Further 
investigation  is  needed  to  determine  the  reason  for  the  bias  in  the 
autocorrelation  function  R(l).   A  general  procedure  for  choosing  the 
proper  signs  to  use  in  the  solution  of  the  nonlinear  algebraic  equations 
is  also  desirable. 

The  Kalman-f ilter  approach  would  be  ai  excellent  method  to  use  in 
the  identification  process  if  a  suitable  sensitivity  function  can  be 
found.   This  function  must  reflect  the  sensitivity  of  the  performance 
index  to  the  changes  in  the  input  distribution  matrix.   A  possible 
first  step  in  finding  the  function  is  to  examine  the  sensitivity  of 
the  measurements  to  the  elements  of  the  matrix. 
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APPENDIX  A 


DERIVATION  OF  THE  AUTOCORRELATION  FUNCTIONS 
FOR  THE  NOISE  ELEMENT 


R(O)    =  EfQ^] 


A 

Vic    '    "2"k-l unwk-n+l    '     vk 


Q     =  d.,w,    +  d„w,     ,   +   ...  +  d..w,       ,  ,  +  v 


R(O)    =«(Vk  +  "2Vl+    •••+dnWk-n-H)2' 

2      2  2  2  2  2  2 

R(O)   =  E  d,  V      +  d/w,     /+...+  d     w.       . ,      +  v. 

Ik  2      k  - 1  nk-n+1  k 

+  2d1d0w.w,     1   +  2d1d„w,  w,     „+...■+■  2d,d  w.  w,        ,  , 
1   2   k  k-1  13k  k-2  Ink  k-n+1 

+  2d2d3wk_lWk_2  +    ...  +  2d2dnwk_lWk_n+1 

+   ...   +  2dlwkvk  +  2d2wk_lVk  +   ...   +  2dnwk_n+1vk] 

»f      2    .   A        2         ,        „      , 
E[wk-i]    =Qw  i¥=°'    l-    •'• 

r.r      2,      A        2 
E[v.  =  a 

k  v 

E[w,     .w.     .]    =  0  i   £  j 

k-i  k-j  J 

,    A 
E[w.     .v.       =  0  i  =  0,    1,    ... 

k-i   k 

R(O)    =  Eld^w/]   +  ■[d2\.12]  +    ...  +  IH,2Vrtll   +  E[vk2) 

R(0)    =  dL2   E[wk2]  +  d/   Etw^]  +   ...   +  E[dn2„k2n+1]  +  a,2 

2     2     2  2      2 

R(0)  =  a   (dn  +  d„  +  ...  +  d  N  +  a 
w   v  1     2  n)v 


(j)  -E[^+jl 


R(j)  -  E[(d  w  +  d  w.  ,  +  ...  +  d  w.      +  v,  )• 
Ik    2  k-1  n  k-n+1    k 

(d,w    +  d  w   '.  +  ...+  d  w.   ,.,-,+  v.,.) 

1  k+j     2  k+]-l  n  k-n+j+1    k+j 
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2  2  2 

R(i)    =  E[d,    w.  w       .    +  d0   w        w  +    .  .  .    +  d      w.        Mw. 

VJ/  Ik  k+j  2      k-1  k+j-1  n     k-n+1  k-n+j+1 

+  d.d.ww  +  d   d   w  w  ~   +    ...    +  d   d   w  w 

12k  k+j-1  13k  k+j-2  Ink  k-n+j+1 

+  d,d„w      ,w  +  d9d„w        w       .        +    ...    +  d   d  w        w.  . 

1   2   k-1  k+j  2    3  k-1  k+j-2  2    n  k-1  k-n+j+1 

+    ...   +  d,d  w       w  +  d0d  w  w         ,.+... 

1   n  k+j   k-n+1  2    n  k+j-1  k-n+1 

+  d      ,d  w.  w.        ,  .   +  v.  ,  .    (d-w.    +  d_w,     .   +    ... 

n-1   n  k-n+1  k-n+j  k+j         Ik  2   k-1 

n  k-n+1  k         1  k+j  2   k+j  - 1  n  k-n+j+1 

+  v   v,  ,  .  ] 
k   k+jJ 

A        2 

E[wk-iwk+j-je]  =  °w  *  =^_j   and   j  =  L>  2j  •••  '  n_1 

E[wk-iwk+j-je]  =  °  i  ^_j   or    j  -  n 

E[Vk+j]  -° 

E[vkwk+j]    =   0  3=0,    1,     ... 

R(l)    =  Eld.d.w.  2  +  d-d-w.  2     +    .  .  .   +  d      _d  w.  2    ,  .  ] 
12k  23  k-1  n-1   n  k-n+1 

R(l)    =  a  2    (d   d     +  dd     +    ...   +  d        d    ) 
w  Li.  i.    3  n-ln 

R(2)   =  E[d.d„w  .      +  d0d.w.  2     +    ...   +  d      9d  w.       .    ] 
13k  2    4  k-1  n-2nk-n+lJ 

R(2)   =  a  2    (d,d0  +  d  d.   +   ...  +  d     _d   ) 
w  1   3  24  n-2    n 


2 
R(n-l)    =   E[dLd2w      ] 

R(n-l)    -   a     d,d 
win 


R(j)    =0  j>  n 


2    n"J 

R(j)   =  a         E     d.d.,  .  j   =   1,    2,    ...,   n-1 

w      .    ,      j    i+i  J 

i=l 


R(j)    -   0  j   >  n 
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APPENDIX  B 


DERIVATION  OF  THE  PREDICTED  MEASUREMENT  VARIANCE 


Vk-rE*l*ti/k-1  =  "*U-i 
<zk-  \/k-i)2  -  'H**k  +  \  -H*ik/k-i'2 
pz  •  E[(zk  -  \/k.i»2i 
pz  =  E[iH*^k-ik/k-i)  +  \i2] 

pz  =  ■[yi*o%[4A.1)  +  vk]    [  (2k  -  ^.p1  H*T  +  vk]  ] 

PZ  =  E[H*   (£    -  ^^    (Zk   -  ^/k_1)T  H*T  +  vk2 

+  H*    (ik   -  ^.l}    vk  +  vk    (ik    -  $k/k_1)T  H*T] 

2    A 
E[vkZ]=  R 

A 
E[vfc]    =  0 

E[(^k  "  ^k/k-i}  (\  '  ^k/k-i)T]  =  pk  (in  steady-state> 

PZ  =  H*  ?v   H*T  +  R 
K 
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appendix  c 

derivations  of  the  sensitivity  functions 

Spz 

The  first  sensitivity  function  to  be  derived  is  -^prr  during  steady- 
state  operation  of  the  filter.   This  is  done  by  finding  an  expression 
for  the  prediction  covariance  during  steady-state  operation,  equating 
elements  on  both  sides  of  the  equation,  and  then  taking  the  desired 
partial  derivatives  of  the  implicit  functions.   The  derivation  follows 

Gk  =  pk/k-iH*T  [H*pk/k-iH*T  +  R1"1 


pk/k  ■  (I  "  GkH*>  pk/k-i 

Pk/k  "   (I   "   Pk/k-l  H*T    [H*  Pk/k-l  H*T  +  R1_1   H*)    Vk'x 

T 
P  ,      =    $*   P    /       $*     +  0 

k+l/k  k/k  v 

r  t  t        -i 

P         ,      =  $*      CI    -    P    ,  H*      TH*   P    ,  H*     +  Rl        H*)    P    / 

k+l/k  LU  k/k-1  L  k/k-1  J  '     k/k-1. 


T 
$*  +Q 


DURING  STEADY-STATE  OPERATION 


P      =  P      =  P 
k+l/k    k/k-1    K 


P  =  $* 
K 


(I  -  P  H*T  [H*P  H*T  +  R]"1  H*)  P. 


K 

FOR  A  SECOND-ORDER  SYSTEM 


KJ 


T 
>*     +  Q 


0      1 

PK  = 

cpx  cp2 

1   0 

3    1 
_    -i 


P       P 

11      12 

P        P 
21      22 


+  R 


[io] 


Q1I  Q12 
Q21  Q22 


P  P 

11  12 

P  P 

21  22 


—                   ™" 

"1 

—         ~ 

Pll   P12 

1 

0     cp1 

P        P 
21     22 

1 

1     cp2 

—           _ 

-> 

—        — 
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PK  = 


-P2i2                                           I      -p 
F^Tl  +  P22  +  Qll            I  P-^TR   Vll^V  +    (cplP21-kP2P22)  +  Q12 
I 

"P21        ,     .  _     .   I  ,      _  .      .     ,  ^^ll       V21, 


(cPtP.,  +cp9P91)  I  (cp.P^-kp^P.,)    (cp. 


P,,+R    ^1    11       ^2   21'      V4Vll    '2    12'    VH"1      P17+R        P..+R' 
11  1111 

+   (cp2P21  +  cp2P22)  +  Q21  i+  cp2    (cp1P21-kp2P22)  +  Q22 


P       =  P 
r21  12 

Q21  =  P12 
EQUATING  ELEMENTS 

"P212 
Pll  =  P— T^  +  P22  +  Qn 

"P21 
P21  =   P12   =  ?—n   (CD1P11  +  ^P  +   (rD2P21  +  CP2P22)  +  Q12 

cp1Pn  cp2P21 

P22    =    (CD1PU   +  CD2P21)     (CP1    -   — —    -    pJYTR> 

+  cp9    (cp,P91  +  cp9P99)  +  Q 


2    VH1   21       ^2   22'        x22 

REARRANGING 

(F.)      Pu2  +  Pn    (R   -    P22    -   N7;L2)   +    (P212    -    P22R   -    NR7l2)    =  0 

2  2 

(G.)     cp2      P21     +  P21    (Pu  +  R   -  cp^)    -    (cp2    P22    Pu  +  cp2RP22 

+   Pu   N7;l72  +  NR7;L72)    =   0 
(H.)      P22    (Pu  +  R   -  cp2      Pn    -  cp2   R)    -    (CP1      PUR   -  CD2    P21 


2  2 

+  2cp1  cp2    P21R  +   PnN72     +  NR72    )    =   0 
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N  =  E    [w      ] 

K. 


Q   =  N 


7i      7i7a 


7172      72 


r*  = 


PZ    =  H*P  H*     +  R 
K 


PZ    =    [1    0] 


PZ   =   Pu  +  R 


P        P 

11      12 


P        P 
r21      22 


+  R 


(F.)      PZ2    -    PZ    (R  +   P22  +  N7L2)   +  P21      =0 

(G.)      cp22    P212  +   P21    (PZ    -  q^R)    -  (cp2P22    PZ  +   PZ    N7l72)    -   0 

2  2  2  2 

(H.)      P22   PZ    (1    -  cp2    )    -  cp1      PZ   R  -    2  cd^2   P21  R  +  cp2      P21 

2      2  2 

+  cp        R      -  PZ   N7        =0 


^PZ 
%7l 


ar 

B7l 


aF 

apz 

^G 
^PZ 

aH 
apz 


aF 
ap21 

aF 
ap22 

ac 
ap21 

aG 
ap22 

aH 
ap21 

aH 
a?22 

aF 
ap21 

aF 
ap22 

ac 
ap21 

ac 
ap22 

aH 
ap21 

aH 
a?22 

J 
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dPZ 


-2N7    PZ 


•N72PZ 
0 


2P 


21 


2r°2      P21  +   PZ    '   ^lR 


■2ro  Q72R  +  2m2      p 


-PZ 

-co2    PZ 

2 

PZ     -    QB2 

PZ 

2PZ    -   R   -    P22    -   Ny 


P21   -  ca2P22    -   N7lr2 

2  2   2  2 

P22-ro2   P22~cpl  R   ~N^2 


2P 


21 


2cn2      P21  +   PZ    "  ®lR 


■2cd  cp2R  +  2ro2      P 


-PZ 

-ro2    PZ 

PZ    -  cd22 

PZ 

^7, 


£f 


ap 


21 


aF 


ap 


22 


ac, 
a7c 


ac 


^p 


21 


ac 

^2  I 


apz 

^7o 


aH 


ap 


21 


aH 

3p 


22 


2P 


21 


■PZ 


apz 

a70 


■N71PZ 
■2N72PZ 


2CD2      P2      +    PZ    -  cp  R 

2 
-2co..cp2R  +  2cd2    P?1 


-CD2PZ 
PZ    -    CD2    PZ 


apz 

ar* 


apz 
apz 

^7o 
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Another  sensitivity  function  of  the  form 

f  (j,  r*)  =  h*  . 


K 


A      .21    <^PZ 

\  Q2i  "  zi/k-i;  J  "  ar* 

i=K 


is  to  be  considered  now.   It  may  be  found  by  expressing  the  filter  gain 
in  steady  state  as  a  function  of  the  steady-state  prediction  covariance. 
The  required  partial  derivatives  may  then  be  found  as  above.   The  follow- 
ing derivation  is  done  assuming  steady-state  operation  of  the  filter. 
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where  -,p .  is  as  given  above 
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